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The dynamics of neutrally buoyant particles transported by a turbulent flow is inves- 
tigated for spherical particles with radii of the order of the Kolmogorov dissipative 
scale or larger The p seudo-penalisation spectral method that has been proposed by 
Pasauetti et all ^00% is adapted to integrate numerically the simultaneous dynamics of 
the particle and of the fluid. Such a method gives a unique handle on the limit of valid- 
ity of point-particle approximations, which are generally used in applicative situations. 
Analytical predictions based on such models are compared to result of very well resolved 
direct numerical simulations. Evidence is obtained that Faxen corrections give dominant 
finite-size corrections to velocity and acceleration fluctuations for particle diameters up 
to four times the Kolmogorov scale. The dynamics of particles with larger diameters is 
dominated by inertial-range physics, and is consistent with predictions obtained from 
dimensional analysis. 



1. Introduction 

A large number of natural and engineering situations involves the transport of spherical 
finite-size particles by a fully developed turbulent flow. This includes the formation of 
planets in the early solar system, rain formation in clouds, the coexistence between several 
species of plankton, and many industrial settings encountered in chemistry and material 
processes. An important feature of such particles is that they do not follow exactly the 
fluid motion but have inertia, a p roperty that leads to the development of inhomogeneities 
in th eir spatial distribution (see Squires fc Eaton 1991 ; iBalkovskv" et a/."200l";'Bec et al\ 
2007t) or to the enhan ce ment of the rate at which they co llide (see Falkovich et al. 2002j; 
Wilkinson e7aII[2006tlZaichik a/.ll2006HBec et aLll2009l ). The modeling of such particles 
generally assumes that their diameter dp is much smaller than the smallest active length- 
scale of the flow, that is the Kol mogorov scale 77, so that they can be approximated by 



points (see Maxev fc Rilev Il983l : ICatignoll Il983 ). Modeling situatio ns where d 



relies on the use of various empirical laws (as reviewed, for instance, in lClift et 



> 



1978h . 



Generally such laws are obtained by considering a particle suspended in a mean laminar 
flow and interacting only through the turbulent wake that it creates, but not with a fully 
developed turbulent environment maintained by an external energy input. 

Recent experimental developments have triggered a renewal of interest in the under- 
stand ing and quantification of finite-size eff e cts in the motion of particles in a turbulent 
flow (|Qureshi et al\ [iooi IVolk et al\ boOSt IXu fc Bodenschat j [20oi ICalzavarini et al 



200i). These works addressed in particular the problem of delimiting the domain of 
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validity of the point-particle models that are largely used in applicative fields, and to un- 
derstand which corrective terms give dominant corrections. Such questions remain largely 
open because of the difficulty in constructing analytically the fluid flow perturbed by the 
presence of the spherical particle. In the following, we briefly review the equations that 
govern the coupled dynamics of the flow and the particle. In an incompressible fluid with 
kinematic viscosity v and subject to an external volumic forcing strain tensor F, the 
velocity field u solves the Navier-Stokes equation 

afM+(M-V)M = -— Vp + i/V^M + V-F, V-u = 0, (1.1) 

which is supplemented by a non-slip boundary condition 

M(Xp(<) + (dp/2)n,i) = l^(i) + ^rjp(t)xn, Vn : |n| = 1 (1.2) 

at the surface dB of the spherical particle. X-p{t) denotes here the trajectory of the center 
of the particle, lp(i) its translational velocity, and r2p(t) its rotation rate. The motion 
of the particle is determined by Newton's second law 

ITT p r 

mp — ^ = (TOp - mf ) g + / V • T dV = (wp - mf ) g + / T ■ AS (1.3) 

Jb Job 

where nip = (7r/6)ppdp is the particle mass (with pp its mass density), mi — (7r/6)pf dp 
the mass of the displaced fluid, g is the acceleration of gravity T — —pl^ + (/x/2) (Vtt + 
'VvJ) + pf¥ denotes the fluid stress tensor, I3 is the identity, and p — pfi/ the dynamic 
viscosity. In addition, the sphere rotation rate flp changes according to the conservation 
of angular momentum 

where n denotes the outward pointing unit-vector normal to the surface and where 
we have assumed that the mass moment of inertia tensor X is that of a uniform solid 
sphere. Solving the system (|l.ip - p.4p is a difficult task as it involves a nonlinear partial 
differential equation for the fluid, which is coupled to a moving boundary condition on 
the sphere. Analytical treatments of such a complex dynamics has only been done when 
neglecting nonlinearities in the flow motion at the scale of the particle, so that (|l.ip 
reduces to the Stokes equation (this leads to the usual point-particle models). 

This study focuses on neutrally buoyant particles, i.e. Pp = pf. This case is of interest 
for applications to problems of plankton dynamics in the ocean or of some types of 
ice crystals in clouds. The goal is here to give a complete description of the dynamical 
properties of particles with sizes of the order of the Kolmogorov dissipative scale rj. The 
paper is organised as follows. In fj2l we consider the model given by the point particle 
approximation. We show that in the case of neutrally buoyant particles, flrst-order flnite- 
size effects are not due to particle inertia but purely stem from Faxen corrections. They 
intervene in the particle dynamics as (dp/A)^, where A designates the Taylor micro-scale. 
These results are validated numerically in fJHlthanks to the use of a new dynamical pseudo- 
penalisation technique that has the advantage of allowing one to use a spectral code in 
order to integrate the Navier-Stokes equation with the proper boundary conditions. 
We show that, both for velocity and acceleration statistics, finite-size effects become 
noticeable for dp ~ 3 77 and that first-order Faxen corrections are relevant up to dp = 
4 77. For dp > 4 r/, the particle dynamics is dominated by inertial-range physics. We 
also present results on acceleration time correlation that confirm this fact. Finally, 21 is 
dedicated to concluding remarks and prospectives. 
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2. Point-particle approximation 

The derivation of point-particle models relies on the assumption that the per turbation 



of the surrounding flow by particles is well des cribed by the Stokes equation (see lGatignol 
1983t iMaxey fc Riievlll983HAuton et a/.lll988l ) . This assumption clearly requires that the 



particle Reynolds number defined with the velocity difference between the fluid and the 
particle is very small. The motion of the neutrally buoyant particle is then given by 

where and account for Faxen corrections. They are averages of the fluid velocity 
over the surface and of the fluid acceleration over the volume of the particle, respectively: 

^■"W^^j/ u{x,t)dS and A'^it)^-^ f ^{x,t)dV, (2.2) 

where D/Dt = dt + u ■ W denotes the material derivative along fluid tracer trajectories. 
The various forcing terms in (|2.ip are, in order of appearance, the combination of the 
inertia force exerted by the undisturbed flow and the added mass, the Stokes viscous 
drag, and the Basset-Boussinesq history force. In the limit when the particle size is 
much smaller than the Taylor micro-scale A, a Taylor expansion of the fluid velocity in 
the vicinity of the particle center leads to 



U^it) = «(Xp, t) + ^dl V2«(Xp, t) + 0[(dp/A)4], 

A^(t) = (DM/Dt)(Xp,<) + ^dl (DV2«/Dt)(Xp,t) + 0[{dp/XY 



(2.3) 



In the limit of particle diameters much smaller than the Kolmogorov scale rj, the Basset- 
Boussinesq history force gives a contribution much smaller than the viscous drag and 
can thus be neglected to leading order. Hence, finite-size neutrally buoyant particles 
obey asymptotically the minimal model equation 

Note that in this model, the particle size enters only the coefficient of the drag force. 
However, as it is now shown, such an effect is actually not sufficient to account for 



leadin g-order corrections due to the particle finite size. Indeed, following iBabiano et 



( 200Cl[ ) and introducing the velocity difference between the particle and the fiuid W{t) 



\^{t) — u{Xp{t),t), one can easily check that 

dW 12;y 

_ = -Ty . Vw(Xp(t), t) - ^W. (2.5) 

This implies that W{t) = eyip{-12v t/dl) rexp[- J* Vu{Xp{s),s) ds] W{0), where Texp 
denotes the time-ordered exponential. Hence the amplitude of the velocity difference 
grows exponentially at large time, i.e. \W{t)\ ~ |VF(0)| exp[— (12i^/(ip -t- Aa)^], where 
A3 is the smallest Lyapunov exponent associated to Texp[— VM(-Xp(s), s) ds]. Because 
of the fiuid fiow incompressibility implying a vanishing sum of the Lyapunov exponents, 
A3 ^ 0. However, when dp is sufficiently small, i.e. when dp ^ •y/l2i7/|A3|, the exponential 
growth rate of \W\ is negative and the particle velocity relaxes to that of the fluid. 
Estimating the value of A3 requires in principle to evaluate the Lyapunov exponents 
associated to the fluid flow strain along particle trajectories. However, because of the 
exponential relaxation of the particle velocities to that of the fluid, these exponents 
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are exactly those computed along tracer tra.iectories. Th e latters have been evaluated 



in direct numerical simulations (see, e.g., Bee et al. 2006f l and, once normalised by the 




inverse of the Kolmogorov eddy turnover time r^, depend weakly upon the Reynolds 
number of the flow: for Re\ varying from 65 to 185, one observes r^Aa = 0.190 ± 2%. To 
summarise, this shows that the minimal model (12. 4p exactly sticks to tracer dynamics 
for sizes smaller than a fixed threshold, i.e. for 



(2.6) 



Hence, for $ ^ $* finite-size effects are not related to inertia but can only stem from 
terms that were neglected in this model. This observation explain why neutrally buoyant 
particles whose dynamics is approximated by the point-particle mo del ()2.4p do not display 
any clustering properties, as observed by ICalzavarini et "01. ( 2008h . 

Let us now estimate the contribution from terms that were neglected, namely from 
the Faxen c orrections and the Basset-Boussinesq history term. For this, we follow the 
approach of Maxev ( 19871 ) developed for small Stokes numbers and write the following 
perturbative Ansatz Vp{t) = u{Xp{t),t) + (dp/A)"/(t) o[((ip/A)"], where the order 
a and the function / are to be determined. Inserting this form in (|2.ip and (|2.3p . one 
obtains that the first-order terms originate from Faxen corrections to the viscous drag, 
so that a — 2 and 



%it)=uiXp{t),t) 



o[id,/xY 



(2.7) 



Note that the synthetic velocity field defined above is divergence-free, implying that no 
effect of particle preferential concentration can be detected with first-order corrections. 
This asymptotic form (|2.7p implies that the particle velocity variance satisfies 



<|VpP)-{M^) 



20 \ 



20 ly 



20 



A 



(2.8) 



where e is the mean turbulent rate of kinetic energy dissipation. As for the variance of 
particle acceleration, one obtains from the time derivative of (|2.7p 



( 


AVp 


■)-{ 


Dm 




DVm 




dt 




Dt 




Dt 



V^Vm 



(2.9) 



where || • || denotes the tensorial Frobenius norm, i.e. ||M|p = trace (M''' M) . Hence we 
expect at small diameters that finite-size effects materialise as a falloff of the two above- 
mentioned dynamical properties of particles that behave quadratically as a function of 
the particle diameter with a coefficient given by Eulerian averaged quantities. 



3. Numerical results using a pseudo-penalisation method 

In order to assess numerically the effect of the particles' finite size on their dynamics, 
a pseudo-spectral code has been adapted to non-trivial geometries by using a pseudo- 
penalisation strategy. The purely spectral part of th e parallel code LaTu, which was 
already used to investigate Lagrangian turbulence by Homann et al\ (2007), solves ac- 
curately the incompressible Navier-Stokes equations. This method treats the evolution 
of the fluid velocity fleld in Fourier space and computes convolutions arising from the 
non-linear terms in physical space. The Fourier transformations are performed by the 
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Rex ""rms £ V 5x St Tj Tr, L Tl 

32 0.17 4.5-10-^ 3-10-^ 1.23-10-2 8 • 10"^ 5 ■ 10"^ 0.8 1.2 6.5 512^ 

Table 1. Parameters of the numerical simulations. Re\ — \J\^m^^^^^Llv-. Taylor- Reynolds num- 
ber, Urms: root-mean-square velocity, e: mean kinetic energy dissipation rate, v: kinematic vis- 
cosity, Sx: grid-spacing, St: time step, r\ = (z/^/e)^''*: Kolmogorov dissipation length scale, 
Tr, = (i^/e)^''^: Kolmogorov time scale, L = (2/3E)'^^^ /£: integral scale, Tl = L/urms: large-eddy 
turnover time, A*'^: number of collocation points. 



PSDFFT-libraryQ The domain is a triple-periodic cube. This method allows for high 
accuracy, precise control of the physical parameters and numerical efficiency. In order to 
maintain a statistically steady state we force the flow by prescribing the energy content 
of the Fourier vectors with moduli 1 and 2. The energy content of each of these two 
shells is kept constant while the individual amplitudes and phases are evolved piecewise 
linearly in time between several random configurations separated by a time 10 T^. The 
advantages of such a forcing are two-fold: it allows one to achieve a statistically isotropic 
large-scale flow and limits the fluctuations to only approximatively 10% of the mean. The 
turbulent characteristics of the flow generated in this way are summarised in Tab.[TJ 

In order to im pose the no-slip bound ary condition on the surface of the spherical 
particle, we follow iPasauetti et M ()2008l) and use a pseudo-penalisation method, which 
consists in imposing a strong drag to the fluid velocity at the particle location, so that 
it relaxes quickly to the particle solid motion. The hydrodynamical forces acting on the 
particle are computed by a Riemann approximation of the integrals appearing in (|1.3p 
and (|1.4p on a homogeneous grid of discrete points located on the surface of the sphere. 
The value of pressure at these points is computed by tri-cubic interpolation. The surface 
integral of the fluid velocity gradient is computed from evaluating the average velocity 
on spherical shells surrounding the particle. At the moment, the simulations are limited 
to a single particle in order to prevent the individual dynamical properties from being 
contaminated by particle-particle hydrodynamical interactions. Notice however that the 
code is very well adapted to situations involving several particles. Because only a single 
isolated particle is considered in the flow, the statistical convergence of particle-related 
quantities requires to perform averages over very long times. Each simulation for a single 
value of the particle diameter required to integrate the flow over more than three hundreds 
large-eddy turnover times. Eight different particle diameters dp are considered within 
the range 2 77 to 14 77. As the pseudo-penalisation technique requires several grid-points 
inside the object in order to correctly impose the boundary conditions, the Kolmogorov 
dissipative scale 77 is resolved with four grid points. This requires the use of double floating 
point precision. Because of the large spatial resolution and the long time integration 
which are required, the simulations are very computationally demanding: this work took 
approximatively four millions of single processor CPU hours. Figure [T] (Left) shows the 
typical vorticity fleld in a cut-plane passing through the center of the particle. Note that 
the signature of a turbulent wake is visible on the right-hand side of the particle. 

To benchmark these simulations, a run with the same parameters as those shown in 
Tab. [1] but without any finite-size particle has been performed. In this simulation, we 
have integrated the motion of passive point particles with a dynamics obeying (|2.4p 
and of passive tracers. It is worthwhile mentioning here that turbulent fluid statistical 
quantities are observed not to depend on the presence of a particle, up to the statistical 



t Parallel 3D Fast Fourier Transforms (P3DFFT) , |http://www.sdsc.edu/us/resources/p3dfft| 
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Figure 1. ie/f; modulus of the vorticity in a slice of the domain that is passing through the 
center of the embedded particle (dark = high vorticity, light = low vorticity); the particle di- 
ameter is here dp = 877. Right: normalised probability density function of the particle velocity 
for various particles sizes, as labeled; the bold dash line corresponds to a Gaussian distribution. 
Inset: deviation of the particle velocity variance ~ (1^1^) from the fluid value, as a func- 
tion of the non-dimensionalised particle diameter $ = dp/r?; the two dashed line represent the 
deviation (|2.8p from the fluid root-mean square velocity that is expected to stem from Faxen 
corrections and a behaviour cx dp''"^, respectively . 



fluctuations due to finite time averages. The particle- free simulation serves thus as a 
reference to estimate tracer statistics and Eulerian averages. 

Figure [T] (Right) shows the probability density function (PDF) of the particle velocity 
components for various particle diameters. Once normalised by their standard deviations, 
these PDFs almost collapse on top of each other and deviate very weakly from a Gaussian 
distribution. The measured variance of the particle velocity that is represented in the 
inset, decreases as a function of the particle size. For small diameters, i.e. for 4> = 
dp/77 < 4, the behaviour of the particle velocity variance is very well described by the 
prediction (|2.8[) obtained from Faxen corrections. For > 4, the deviation from the 
fluid velocity variance behaves as ^-^^^ . This power-law is dimensionally compatible with 
Kolmogorov 1941 scaling and indicates that particles with such diameters respond to the 
inertial-range physics of turbulence. Noticeably, in the low-Reynolds-number flow that 
we are considering here, there is no inertial range in the sense usually defined through 
velocity scaling properties. Hence it seems that particle dynamical properties are much 
more amenable to dimensional estimates than fluid turbulent quantities. 

We next turn to particle acceleration statistics. Figure [2] (Left) represents the compo- 
nentwise normalised variance of the particle acceleration oq = {{dVp /dty)e~'^^^r]'^^^ as 
a function of the non-dimensionalised diameter $ = dp/rj, both for the real particles as 
well as for the minimal point model ()2.4p . For real spherical particles, one distinguishes, 
as in the case of velocity variance, between two behaviours. When <I> = dp/rj < 4, finite- 
size effects in the acceleration variance are very well captured by Faxen corrections and 
are very close to the prediction (|2.9p . Note that, thanks to their isotropic form, the 
sub-leading terms appearing in (|2.9p could be evaluated here through the pressure and 
velocity spect ra. When $ > 4, an inertial-range behaviour with oq cx dp ''^^ is attained. 



As argued by lQureshi et al. ( 2007 ). the variance of finite-size particle acceleration is re- 
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Figure 2. Left: acceleration variance of real spheres and of point particles as a function of 
the non-dimensionalised particle diameter <1> = dp/77; error bars correspond to an estimation of 
the standard deviations. The two dotted vertical line indicate the Taylor micro-scale A and the 
critical value (|2.6|) above which point particles deviate from fluid tracers. The dashed curve cor- 
responds to the behaviour (|2.9p predicted from Faxen corrections. Right: normalised probability 
density function of the component-wise particle acceleration for various particle diameter, as 
labeled. The bold dashed line on whic h data almost collaps e corresponds to the log-normal fit 
l|3.1|) of experimental data proposed bv lQureshi et al\ (|2007l '). 



lated to that of the fluid pressure integrated over a sphere of diameter dp. The power 
—4/3 that is observed here differs from the value —2/3, which was measured by Qureshi 
et al. However, as already stressed for instance in Bee et a/] (|2007.) , pressure scaling in 
low Reynolds number flows is often dominated by sweeping, leading to a behaviour of 
the pressure increments \p{x -\- () — p{x)\ ~ f-^^ . While no scaling of pressure can be 
detected in the present simulations, the robustly observed —4/3 law can be explained 
with such a sweeping-dominated pressure spectrum. Another observation is that numer- 
ics confirm the presence of the threshold (|2.6p predicted in the previous section for the 
minimal point particle model: indeed the numerical integration of point particles obeying 
(|2.4p shows that when <i> < 8, the acceleration variance of the latters is undistinguish- 
able from that of tracers. For $ > 8, the point particle model gives an enhancement of 
acceleration, which is incompatible with measurements done with real particles at the 
Reynolds number considered here. This stresses the irrelevance of such a model in the 

case of neutrally buoyant particles. Note finally that in the case of tracers, the constant 

1/2 I 

an is known to show a Reynolds number dependence oq oc (see, e.g., iVoth et al. 

2002I ). The measured value of approx imatively 1.3 i s in go od agreement with the value 
that was measured experimentally by Qureshi et al. ( 2007 ). 

Figure [5] (Right) represents the PDF of acceleration components norm alised to unity 
varian ces for various values of the particle diameter. As already stressed in I Qureshi et al 



(12007|), the dependence upon $ = rfp/77 is very weak. Data can be fitted by the function 



p{a)= [exp(3sV2)/(4\/3)] {l-erf [(ln|x/\/3|+2s2)/(\/2s))}, (3.1) 



which was proposed bv lMordant et al. (2004) for the acceleration PDF of fluid tracers. 
Numerical results almos t collapse to such a distribution with a value of the parameter 
s = 0.62, as observed bv lQureshi et al. ( 2007 ). 

Other measurements relate to two-time statistical properties of particles. Figure [3] 
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Figure 3. Component- wise acceleration time correlation C(r) defined in (|3.2|l for various par- 
ticle sizes, as labeled. Inset: integral correlation time Ti defined from (|3.3p as a function of the 
non-dimensionalised particle diameter $ = d^/ri. 



represents the acceleration time correlation 

as a function of the time lag r for various values of the particle diamet er. The numerical 



measu rements reported here are in good agreement with the results of ICalzavarini et 



( 20091 ) ■ Surprisingly one observes that C(t) deviates only very weakly from the tracer ac- 



celeration temporal correlation for diameters less than A-q, that is when Faxen corrections 
are expected to be of relevance to capture first-order finite-size effects. This behaviour is 
even clearer when looking at the diamet er dependence o f the correlation time for particle 
acceleration. For this we follow Calzav arini et al. ( 2009f) and introduce the integral time 



Ti= r C(t) dr, (3.3) 

JO 

where Tq is the first zero-crossing time. The inset of Fig. [3] represents Ti/r,, as a function 
of $ = d-p/r}. When $ < 4, this integral correlation time is, up to numerical errors, 
undistinguishable from the value obtained for tracers. When > 4, the correlation 
time increases much faster as a function of $ and follows approximatively the power-law 
behaviour Ti ~ $2/3^ This indicates again that turbulent inertial physics is pulling strings 
at such values of the particle diameter, and that the relevant time scale is then given by 
the eddy turnover time e~^^'^d^^^ associated to the particle size. 



4. Concluding remarks 

In this paper it is shown that first-order finite-size corrections to the dynamics of neu- 
trally buoyant particles are due to Faxen terms and not to particle inertia. This leads to 
several predictions on turbulent velocity and acceleration second-order statistics. Using a 
pseudo-penalisation spectral method, these predictions have been confirmed numerically 
for particles with diameters dp up to 4 rj. Higher-order statistics of velocity and accel- 
eration seem much less sensitive to the finiteness of the particle sizes, since numerical 
observation suggest that, once normalised to unit variance, their PDFs collapse on top 
of each other for various values of the particle diameter. 
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The irrelevance of particle inertia with respect to Faxen terms at small particle diame- 
ters has noticeable consequences. First, it implies that the particle dynamics is very well 
approximated by the advection by a synthetic flow, which is incompressible to both the 
leading and the first orders. This means that the effect of preferential concentration onto 
the dynamics of neutrally buoyant particles is very weak. Secondly, such an observation 
clearly questions the relevance of inertial-particle models for density contrasts between 
the particle and the fluid that are close to one. A third consequence is related to the 
fact that corrective terms apply when the particle diameter is much smaller than the 
Taylor micro-scale A rather than the Kolmogorov dissipative scale 77. This fact might 
partly explain the difficulties in matching experiments and numerical model including 
Faxen corrections (see, e.g.. ,Calzavarini et al. 2009.). 

For particle diameters larger than Arj, inertial-range physics comes into play. A 
striking observation is the relevance of the dimensional estimates that are given by Kol- 
mogorov 1941 theory, even when the fiuid flow Reynolds number is so low that no scal- 
ing range can be observed for Eulerian velocity statistics. Understanding the reasons of 
such a behaviour requires to investigate in a more systematic manner the dynamics of 
particles with inertial-range sizes in fully developed turbulent flow. Applying the pseudo- 
penalisation method to larger particles and higher fluid flow Reynolds numbers is the 
subject of on-going work that is mainly focusing on describing the flow modification 
induced by the presence of the spherical particle. 
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